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We consider insulating phases of cold spin-1 bosonic particles with antiferromagnetic interactions, 
such as 23 N a, in optical lattices. We show that spin exchange interactions give rise to several distinct 
^ ■ phases, which differ in their spin correlations. In two and three dimensional lattices, insulating 

' phases with an odd number of particles per site are always nematic. For insulating states with an 

even number of particles per site, there is always a spin singlet phase, and there may also be a first 
order transition into the nematic phase. The nematic phase breaks spin rotational symmetry but 
preserves time reversal symmetry, and has gapless spin wave excitations. The spin singlet phase does 
not break spin symmetry and has a gap to all excitations. In one dimensional lattices, insulating 
^ — . phases with an odd number of particles per site always have a regime where translational symmetry 

is broken and the ground state is dimerized. We discuss signatures of various phases in Bragg 
scattering and time of flight measurements. 
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O ' I. INTRODUCTION 

•4— » . 

Modern studies of quantum magnetism in condensed matter physics go beyond explaining details of particular 
experiments on the cuprate superconductors, the heavy fermion materials, organic conductors, or related materials, 
and aim to develop general paradigms for understanding complex orders in strongly interacting many body systems 
pi I2I 13. Irl I7L I8ll3. ITpj . Spinor atoms in optical lattices provide a novel realization of quantum magnetic systems 
that have several advantages compared to their condensed matter counterparts, including precise knowledge of the 
underlying microscopic models, the possibility to control parameters of the effective lattice Hamiltonians, and the 
absence of disorder. 

Degenerate alkali atoms are generally considered as a weakly interacting gas due to the smallness of the scattering 
length compared to the inter particle separation The situation may change dramatically either when atomic 

scattering length is changed by means of Feshbach resonance |T3 | , or when an optical potential created by standing laser 
beams confines particles in the minima of the periodic potential and strongly enhances the effects of interactions. In the 
£S) \ latter case existence of the nontrivial Mott insulating state of atoms in optical lattices, separated from the superfluid 
*vO • phase by the quantum phase transition (SI transition), was demonstrated recently in experiments |l3l Il4l Hjj. Low 
energy (temperature) properties of spinless bosonic atoms in a periodic optical potential are well described by the 
CO , Bose-Hubbard Hamiltonian [la| 

O . 

Hbh = -t ^(ajoj + ptqj) - ij, fij + -y hi(hi - 1), (1) 

£h (ij) 1 i 

Parameters of Q may be controlled by varying the intensity of laser beams, so one can go from the regime in which the 
kinetic energy dominates (weak periodic potential, t >> Uo), to the regime where the interaction energy is the most 
important part of the Hamiltonian (strong periodic potential, t << Uo). For integer fillings (number of atoms per 
lattice site), the two regimes have superfluid and Mott insulating ground states, respectively, as can be obtained from 
the mean-field analysis of the Bose-Hubbard Hamiltonian |0, Il7| . In the superfluid phase, atoms are delocalized in 
the lattice, fluctuations in the number of atoms in each site are strong, and there is a phase coherence between different 
sites. In the insulating state, atoms are localized, fluctuations in the particle number at each site are suppressed, and 
there is a gap to all excitations. Such an insulating state represents a correlated many body state of bosons, where 
strong interactions between atoms result in a new ground state of the system. 

In conventional magnetic traps, spins of atoms are frozen so effectively that they behave like spinless particles. In 
contrast, optically trapped atoms have extra spin degrees of freedom which can exhibit different types of magnetic 
orderings. In particular, alkali atoms have a nuclear spin I = 3/2. Lower energy hyperfine manifold has 3 magnetic 
sublevels and a total moment S — 1 . Various properties of such condensate in a single trap were investigated 120L 
[2ll I22I I23I l24l]. For example, for particles with antiferromagnetic interactions, such as 23 Na, the exact ground state 
of an even number of particles in the absence of a magnetic field is a spin singlet described by a rather complicated 
correlated wave function [2j|. However, when the number of particles in the trap is large, the energy gap separating 
the singlet ground state from the higher energy excited states is extremely small, and for the experiments of ref. 
[l8[ . the precession time of the classical mean-field ground state is of the order of the trap lifetime. So, experimental 
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FIG. 1: General phase diagram for S — 1 bosons in 2d and 3d optical lattice. Detailed discussion of the phase diagram, 
including explicit expressions for various phase boundaries, is given in Section EH 

observation of the quantum spin phenomena in such systems is very difficult. To amplify quantum spin effects one 
would like to have a system with smaller number of particles and stronger interactions between atoms. Hence it is 
natural to consider an idea of S = 1 atoms in an optical lattice, in which one can have a small number of atoms per 
lattice site (in experiments of Ref. |l4j this number was around 1-3) and relatively strong interactions between atoms. 

In this paper we study bosonic 5 = 1 atoms in optical lattices with spin symmetric confining potentials and 
antiferromagnetic interaction between atoms. We demonstrate that spin degrees of freedom result in a rich phase 
diagram by establishing the existence of several distinct insulating phases, which differ from each other by their spin 
correlations. 

In the insulating state of bosons in an optical lattice fluctuations in the particle number on each site are suppressed 
but not frozen out completely. Virtual tunnelling of atoms between neighboring lattice sites gives rise to effective 
spin exchange interactions that determine the spin structure of the insulating states (spin exchange interactions for 
5=1/2 bosons in optical lattices were discussed previously in j27ll28l|V 

We will show that in two and three dimensional lattices insulating states with an odd number of atoms per site are 
always nematic, whereas insulating states states at even fillings are either singlet or spin nematic [25|, depending on 
the parameters of the model. In one dimensional systems even more exotic ground states should be realized, including 
the possibility of a spin singlet dimerized phase that breaks lattice translational symmetry [3ll l32| . The 2d and 3d 
general phase diagram, including singlet, nematic and supcrfluid phases, is shown in Fig. ^ The extended version of 
this diagram, including discussion of various transition lines, is presented in section IVTl 

It is useful to point out that the lattice model for spin-1 bosons, which we analyze here, is very general and may 
also be applicable to systems other than cold atoms in optical lattices. For example, triplet superconductors in strong 
coupling limit may be described by a similar Hamiltonian, and some of the phases discussed in this article may 
correspond to non-BCS states of such superconductors [33| . 

The paper is organized as follows. In section[H]we provide derivation of the Hubbard-type Hamiltonian for spin-1 
bosons in optical lattices starting from microscopic interactions between atoms, and describe some general properties 
of our model. In section lTTTI we derive an effective spin Hamiltonian which is valid for any odd number of atoms per site, 
N, in the limit of small tunnelling between sites. We demonstrate equivalence between our system and a Heisenberg 
model for S = 1 spins on a lattice with biquadratic interactions and argue that the ground state is a nematic in two 
and three dimensions and is a dimerized singlet in Id. In section Hvl we derive effective spin Hamiltonian for a system 
with N = 2 atoms per site, valid deep in the insulating regime, and use mean-field approximation to determine the 
phase boundaries between isotropic and nematic phases. In section [3 we derive effective spin Hamiltonian for the 
limit of large number of particles per site N » 1 and small tunnelling, and discuss isotropic- nematic transition for 
even N. In section WH we summarize our results and review the global phase diagram for spin-1 bosons in optical 
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lattices. Finally, in section IVHI we discuss approaches to experimental detection of singlet and nematic insulating 
phases of S = 1 bosons. Details of technical calculations are presented in Appendices A-D. 

II. DERIVATION OF BOSE-HUBBARD MODEL FOR SPIN-1 PARTICLES 

At low energies scattering between two identical alkali atoms with the hyperfine spins S = 1 is well described by 
the contact potential 01 

V"(rx-ra) = 5^ - r 2 )(g V + g 2 V 2 ), (2) 
55 = 4nh 2 a s /M. (3) 

Here Vs is the projection operator for the pair of atoms into the state with total spin 5 = 0,2; as is the s-wave 
scattering length in the spin S channel; and M is the atomic mass. When writing |(5J we used the fact that s-wave 
scattering of identical bosons in the channel with total spin 1 is not allowed by the symmetry of the wave function. 
Interaction J2J) can be written using spin operators as 

V{tx -r a ) = 6{t\ -r 2 )( + — - — SiS 2 ). (4) 

For example, in the case of 23 Na, g 2 > go, and we find effective antiferromagnetic interaction, as was originally 
discussed in [ljj |2(| ■ 

Kinetic motion of ultracold atoms in the optical lattice is constrained to the lowest Bloch band when temperature 
and interactions are smaller than the band gap (this is the limit that we will consider from now on). Atoms residing 
on the same lattice site have identical orbital wave functions and their spin wave functions must be symmetric. If we 
introduce creation operators, a[ , for states in the lowest Bloch band localized on site i and having spin components 
(j = { — 1,0,1}, we can follow the approach of |T(j and write the effective lattice Hamiltonian 

H=-t^2 ( a L a j* + a )a a i<?) + -J- ^2ni(ni - 1) + -— ' £(S? - 2n t ) - /i^n*, (5) 



<«>■ 



where 



= £4 



(6) 



is the total number of atoms on site i, and 



Si = ^aj .T . (7 /aj -' (7) 



is the total spin on site i (T aa > are the usual spin operators for spin 1 particles). The first term in (JSJ describes spin 
symmetric tunnelling between nearest-neighbor sites, the second term describes Hubbard repulsion between atoms, 
and the third term penalizes non-zero spin configurations on individual lattice sites. The origin of this spin dependent 
term is the difference in scattering lengths for S — and S — 2 channels as was discussed in j2^|. Finally, the fourth 
term in JSJ is the chemical potential that controls the number of particles in the system. 

Hamiltonian JSJ carries important constraints on possible spin states of the system. The first of them derives from 
the fact that the total spin of a system of N spin-1 atoms cannot be bigger than N, so for each lattice site we have 

Si < N t . (8) 

The second constraint is imposed by the symmetry of the spin wave function on each site 

Si + Ni = even. (9) 

Optical lattices produced by far detuned lasers with wavelength Xi — 2ir/\ki\ create an optical potential V(r) — 
'^2 i Vi sin 2 ki r, with fcj being the wave vectors of laser beams. Using various orientations of beams, one can construct 
different geometries of the lattice. For the simple cubic lattice, parameters of (JSJ can be estimated as 

U 2 = 
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FIG. 2: U2 and t for Na atoms in the simple cubic optical lattice created by three perpendicular standing laser beams with 
A = 985nm. Vo is the strength of the optical potential and En = h 2 k 2 /2M is the recoil energy. The ratio of the interaction 
terms in |J5J, U2/U0, is fixed by the ratio of the scattering lengths and is independent of the nature of the lattice (U2/U0 ~ 0.04 
for 23 Na) 
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where Er = h 2 k 2 /2M is the recoil energy and x — Vo/Er. Note that the ratio U2/U0 is fixed by the ratio of 
scattering lengths, 02/00, for all lattice geometries. Scattering lengths for 23 Na given in [2£| are ai = (52±5)a,B and 
ao = (46 ± 5)as, where as is the Bohr radius. This corresponds to < 02 — ao << 2a2 + ao, so the spin dependent 
part of the interaction is much smaller than the spin independent one. Throughout this paper we will always assume 
< U2 « Uo- While applying results of this paper for the case of 23 Na, one should note that errors in the estimation 
of the exact value of U2/U0 are very big. While considering the spin structure of Mott insulating phases, we will 
assume that U2/U0 is small enough to see the interplay between tunnelling and spin dependent U2 term before the 
superfluid-insulator transitions take place. The positions of superfluid-insulator transitions and the validity of this 
assumption will be discussed in detail in section IVTI We will use the value U^/Uq = 0.04 to make estimates of various 
phase boundaries. In Fig. [3 we show U2/T1 and t/h as a function of the strength of the optical potential for a three 
dimensional cubic lattice produced by red detuned lasers with A = 985nm. 

Superfluid-insulator transition is characterized by a change in fluctuations in particle number on individual lattice 
sites. When the spin dependent interaction (U2) is much smaller than the usual Hubbard repulsion (Uo), the superfluid 
- insulator transition is determined mostly by Uo. The spin gap U% term, however, is important inside the insulating 
phase, where it competes with the spin exchange interactions induced by small fluctuations in the particle number, 
and an interesting spin structure of the insulating states appears as a result of such competition. The spin structure 
of the insulating phases of spin-1 bosons in optical lattices will be explored in this paper. 

In what follows we will often find it convenient to use particle creation operators that transform as vectors under 
spin rotations. Such representation may be constructed as 



t _ t t 

a z — a 0' a x 



(a-) f -ai + .(a_)t+a 1 , 

L n ' — 



V2 



V2 



(10) 



Operators a,t x ,y,z} satisfy the usual bosonic commutation relations, and they can be used to construct spin operators 
as 



S:, 



i&abcQ'ibQ'ic-i S^ — (&bji3ym &bm ^771)^5 ^7^^^^ 



(11) 



We can verify the transformation properties of ai x y ,z} by noting that 



[S a ,ab] 



(12) 



[Sa, a\] 
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(13) 
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Using these operators the hopping term in the Hamiltonian (JU may be rewritten as 

~t H ( a lp a 3 P + a ]p a i P ) 

<ij>,p£{x,y,z} 

and it is invariant under global spin rotations. We will use this property later to simplify calculations and classify 
eigenstates of effective interaction by the total spin of two neighboring sites. 



III. INSULATING STATE WITH AN ODD NUMBER OF ATOMS 



A. Effective Spin Hamiltonian for small t 

We start with the insulating state of the Hamiltonian (JjjJ with an odd number (N = 2n + 1) of bosons per site in 
the limit t — 0. The number of particles on each site is fixed, and the bosonic symmetry of the wave function requires 
that the spin in each site is odd. The interaction \J% term is minimized when the spins take the smallest possible 
value 5=1. In this limit the energy of the system does not depend on the spin orientations on different sites. When 
t is finite but small, we expect that we still have spin 5 = 1 in each site, but that boson tunnelling processes induce 
effective interactions between these spins. In this section we will compute such interactions in the lowest (second) 
order in t. We will also discuss conditions for which our effective Hamiltonian provides an adequate description of the 
system. 

In the second order perturbation theory in t, we generate only pairwise interactions between atoms on neighboring 
sites, so we can write the most general spin Hamiltonian for 5 = 1 particles that preserves spin 50(3) symmetry 

H = - Jo - Ji SiSj - J 2 &&i) 2 - ( 14 ) 

(ij) <ij> 

Here (ij) labels near neighbor sites on the lattice. Absence of the higher order terms, such as (SiSj) 3 , follows from 
the fact that the product of any three spin operators for an 5 = 1 particle can be expressed via the lower order terms. 

To find the exchange constants Jo, 1,2, we need to consider virtual processes that create a state with N = 2n, Nj = 
2n + 2, and Ni = 2n + 2, Nj = 2n. The difference in energy between the intermediate state and low energy Si = Sj = 1 
subspace is of order Uq, and while our subspace is much lower in energy, the second order perturbation theory is valid. 

It is convenient to rewrite the Hamiltonian (|14|) as 

U = e ^P ij -(0)+e 1 ^Py(l) + e 2 ^P ii (2), (15) 

<y> (a) ^i) 

eo = -4J 2 + 2Ji-J , 
ei = — J2 + Ji — Jo, 

e 2 = -J2-J1-J0. (16) 

Here Py (S) is a projection operator for a pair of spins on near neighbor sites i and j into a state with total spin 
Si + Sj — S (S = 0, 1, 2). Equivalence of (|14fl and (|15|) can be proven by noting simple operator identities for two 
spin one particles 

1 = Py(0)+P ii (l)+P iJ (2), 

{Si + Sj) 2 = 4 + 25,5, = 2P tf (l) + 6Pij (2), 

(Si + Sj) 4 = 16 + 165,5,- + 4(5i5j) 2 = 4P y - (1) + 36Py(2). (17) 

Note that states |5j = 1, Sj — 1; Si + Sj = 5) have only the trivial degeneracy corresponding to possible projections 
of total spin 5 on a fixed quantization axis D$ = 25 + 1. 

Since we know a general form of our effective Hamiltonian, we can compute £0,1,2 by calculating the expectation 
values of energy for arbitrary states in the appropriate subspaces 

t 2 I < m \(4 P a jP + a]p ai p)\ Si = 1; S j Z 1; Si + Sj S > I 2 (18) 
■~ E m — Eq 
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FIG. 3: Ratio J1/J2 for the effective spin Hamiltonian 11141 for an odd number of bosons N = 2n + 1 per site. U2 — 0.04£/o- 



Here Eq = 2U 2 is the energy of the configuration with N = 2n + 1 bosons in each of the two wells, and E m is the 
energy of the intermediate (virtual) states, m, that have 2n and 2n + 2 bosons in the two wells, respectively. Both 
energies should be computed in the zeroth order in t. 

It is useful to note that the tunnelling Hamiltonian is spin invariant; therefore, intermediate states in m summation 
in i|18fl should also have total spin S. Another constraint on the possible states m comes from the fact that the 
tunnelling term can only change the spin on each site by ±1 since, in a Hilbert space of each well, operators a\ p , di P 
act as vectors, according to their transformational properties 1)12 (1 -1)1,' ^ . 

Direct calculations in Appendix lAl give 

4t 2 (n + l)(2n + 3) 16i 2 n(5 + 2n) 
€ °~ 3(C/ - 2U 2 ) 15(U + 4U 2 ) ' ( } 

_ 4t 2 n(5 + 2n) 

61 -5(U + 4U a y m 

28i 2 n(5 + 2n) 4(15 + 20n + 8n 2 ) 
&2 ~~ 75(J7 + 4£/ 2 ) 15(C/ + U 2 ) ' ( ' 



Combining ^SJ-lEU w ^ find 



Jo 


4(15 + 20n 


+ 8n 2 ) 


4(l + n)(3 + 2n) 


t 2 


45(?7 + 


u 2 ) 


9(U + 2U 2 ) 


Jl 


2(15 + 20?i 


+ 8n 2 ) 


16(5 + 2n)n 


t 2 


15((7 + 


u 2 ) 


75((7 + 4C/ 2 )' 


■h 


2(15 + 20n 


+ 8n 2 ) 


4(1 + n)(3 + 2n) 


t 2 


45(C/ + 


u 2 ) 


+ 9(U - 2U 2 ) 



128(5 + 2n) 
225(U + 4U 2 ) ' 

4n(5 + 2n) 
225(C/ + 4U 2 ) ' 



(22) 



It will turn out that the ratio between J\ and J 2 determines magnetic ground state, and its dependence on n is 
quite fast, as shown in Fig. [3] 

We now discuss limitations of the Hamiltonian l|14l) with l|22l) . In the insulating state with exactly one boson per 
site, near neighbor interactions always have the form 1)14)1. Explicit expressions for the J's given in (1221) only apply in 
the limit t << Uq. When t becomes comparable to Uq (but we are still in the insulating phase), higher order terms 
become important, including the possibility of spin coupling beyond the near neighbor sites. In the insulating state 
with more than one boson per site (JV = 2n + 1, n > 0), we have an additional constraint: we should be able to 
neglect configurations with spins on individual sites higher than 1. Matrix elements for scattering into such states are 
of the order of (Nt) 2 /Uo (see (I22H . and their energy is set by U 2 . Therefore, the Hamiltonian 1)14)1 applies only when 
Nt « (UqU^ 1 / 2 , which is well within the insulating state when U 2 « U Q (SI transition takes place for Nt ~ Uq). 
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Stot 


S1S2 (S1S2) 2 


Energy 





-2 4 


2Ji - 4J 2 


1 


-1 1 


Ji- Ja 


2 


1 1 


-Ji - J2 



TABLE I: Eigenstates of a two site problem I23H 



B. Phase diagram 

To understand the nature of the Hamiltonian (|14f> in the relevant regime of parameters J 2 > J\ > 0, it is useful to 
start by considering a two-site problem 

W12 = -JiS x S 2 - J 2 (S 1 S 2 ) 2 (23) 
with Si — S 2 = 1. Eigenstates of lj2l?|l can be classified according to the value of the total spin Stot, and their 
energies may be computed using 2S\S 2 — Stot{Stot + 1) — 4. Two spin one particles can combine into Stot = 0, 1, 
and 2. The J± term in (|23|1 favors maximizing SiS 2 by making the fully polarized Stot = 2 state. By contrast, the 
J 2 term favors maximizing (SiS^) 2 by forming a singlet state Stot = (see Table QJ. So, the latter term acts as an 
effective antiferromagnetic interaction for this spin one system, and it dominates for J 2 > J\. If we go beyond a two 
site problem and consider a large lattice, we see that each pair of near neighbor sites wants to establish a singlet 
configuration when J 2 > J\. However, because one cannot form singlets on two different bonds that share the same 
site, some interesting spin order, whose precise nature will depend on the lattice and dimensionality, will appear. 



1. Phase diagram for d — 1 

^From the discussion above we see the conflict intrinsic to the Hamiltonian (JT2J: each bond wants to have a singlet 
spin configuration, but singlet states on the neighboring bonds are not allowed. There are two simple ways to resolve 
this conflict: 

A) Construct a state that mixes 5 = and S — 2 on each bond but can be repeated on neighboring bonds; 

B) Break translational symmetry and favor singlets either on every second bond. 
At the mean-field level, solution of the type A is given by 

\N) = Y[\Si = l,mi = 0). (24) 

i 

This can be established by noting that for any neighboring pair of sites we indeed have a superposition of S — and 
5 = 2 states 

\Si = l,m t = 0}\Sj = l, mj = 0) = -i I Stot = 0) + \j\\Stot = 2,m tot = 0). (25) 

State 1|24[) describes a nematic state that has no expectation value of any component of the spin (S*' y ' z — 0), but spin 
symmetry is broken since ((5f) 2 ) = ((5f) 2 ) = 1/2 and ((S*) 2 ) = 0. It is useful to point out the similarity between 
wave function l|25|) that mixes singlet and quintet states on each bond, and a classical antiferromagnetic state for spin 
1/2 particles that mixes spin singlets and triplets on each bond. Coleman's theorem [3(j (the quantum analog of 
Mcrmin- Wagner theorem) forbids the breaking of spin symmetry in d — 1, even at T — 0. However, a spin singlet 
gapless ground state that has a close connection to the nematic state (|24H has been proposed in |39l |41| for J 2 close 

to J\. 

The simplest way to construct a solution of type B is to take 

\D) = IJ \S i = l,S i+1 = l,S i +S i+1 = 0). (26) 

i— 2n 

Such a dimerized solution has exact spin singlets for pairs of sites 2n and 2ra + 1, but pairs of sites 2n and 2n — 1 are 
in a superposition of 5 = 0, 1, and 2 states. 

According to the variational wave functions l|24l) and (|2(j|> , the dimerized solution becomes favorable over a nematic 
one only for J 2 /J\ > 3/2 in d = 1. However, numerical simulations [38| showed that for J 2 > J%, the ground state 
is always dimerized. It is a spin singlet and has a gap to all spin excitations. This means that the variational wave 
function l|26|) may only be taken as a caricature of the true ground state, although it captures such key aspects of it, 
such as broken translational symmetry and the absence of spin symmetry breaking. 
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2. Phase diagram for d = 2, 3 

The nematic state for the Hamiltonian (|14f> in a simple cubic lattice (d = 3) for J 2 > J\ has been discussed using 
mean-field calculations[34|],a semiclassical approach[3^, and numerically (3(j. Finally, recent work of Tanaka et.al j37| 
provided a rigorous proof of the existence of the nematic order at least in some part of this region, which satisfies 
2.66 Ji > J2 > 2 J\. The variational state for the nematic order may again be given by equation (|24|l and its mean field 
energy is E^ F = — 2J 2 - It is important to emphasize, however, that the actual ground state is sufficiently different 
from its mean- field version H24JI . It is possible to write down dimerized states with energy expectation lower than 
— 2J2; however, numerical results j3|| suggest that the ground state doesn't break translational symmetry. The way 
to obtain a more precise ground state wave function is to include quantum fluctuations near the mean field state, as 
was done in |35j . Hence, the mean- field wave function l|24|l does not provide a good approximation of the ground state 
energy of the nematic state. Nevertheless, it is useful for the discussion of order parameter and broken symmetries of 
the nematic state. 

In the nematic state, spin space rotational group 0(3) is broken, though time reversal symmetry is preserved. The 
order parameter for the nematic state is a tensor 

Qab={S a S b )- S -f{S 2 ). (27) 

In the absence of ferromagnetic order (S a Sb) — (5f,5 a ); hence, Q a b is a traceless symmetric matrix. The minimum 
energy of Ijl4|) is achieved for Q a b that has two identical eigenvalues, which corresponds to a uniaxial nematic 
Then, the tensor Q a b can be written using a unit vector d as 

Qab = Q(d a d b - -<U). (28) 

Vector d is defined up to the direction (i.e. ±d are equivalent) and corresponds to the director order parameter 1 121 . 
For the mean-field state (|24|l . the director d can also be defined from the condition that locally our system is an 
eigenstate of the operator dS with eigenvalue zero. However, such a definition may not be applied generally. 

The nematic phase behaves in many aspects as antiferromagnetic^^i the direction of d*being analogous to staggered 
magnetization. Namely in weak magnetic fields, d aligns itself in the plane perpendicular to magnetic field, and spin- 
wave excitations have linear dispersion j^, with velocity 



c= V2zJ 2 (J 2 - Jx). 

Nematic phases for the system of spin-1 particles have been considered before in literature [3^l3^l3^l3^l3^l3^.l40ll4lj 
so we will not discuss them here more extensively. 



IV. INSULATING STATES WITH TWO ATOMS PER SITE 



In this section we consider an insulating state of two bosons per site. Possible spin values for individual sites are 
5 = and S = 2. In the limit t — 0, the interaction part of the Hamiltonian (the U 2 term) is minimized when 5 = 0. 
The amplitude for creating 5 — 2 states, as well as the exchange energy of the latter, is of the order of t 2 /Uq. So, 
when t is of the order of (C/0^2) 1 ^ 2 or larger, we may no longer assume that we only have singlets in individual sites, 
and we need to include 5 — 2 configurations in our discussion. This regime is still inside the insulating phase for small 
enough U2/U0 (the superfluid-insulator transition takes place for zt ~ Uq). In this section we will assume that U 2 /Uq 
is small enough, so that 5 = 2 becomes important in the insulating phase, before the transition to superfluid. More 
careful consideration of superfluid transition line and comparison with the case of 23 Na will be presented in section 
IVII In section llV Al we exactly solve the problem for two wells. In section llV Bl we derive an effective Hamiltonian 
that takes into account competition between spin gap of individual sites, that favors 5 = everywhere, and exchange 
interactions between neighboring sites that favor proliferation of 5 = 2 states. Mean field solution of the effective 
magnetic Hamiltonian is considered in section TlV Cl and we find first order quantum phase transition from isotropic to 
nematic phase. We discuss collective excitations in sections llV Dl and IIV El and the effects of magnetic field in section 
IIV Fl We note that the state with N — 2 has an advantage over states with higher N from an experimental point of 
view since it has no three-body decays. 
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A. Two site problem: exact solution 



To construct an effective magnetic Hamiltonian for this system, we note that in the second order in t it can be 
written as a sum of interaction terms for all near neighbor sites (identical for all pairs of sites). These pairwise 
interactions can be found by solving a two well problem and finding the appropriate eigenvalues and eigenvectors in 
the second order in t. 

The Hilbert space for two sites with two atoms in each well is given by the direct sum of the following subspaces: 

S 2 = 0,Si + S 2 = 0>, 
S 2 = 2,S 1 +S 2 = 0>, 
0,S 2 = 2,S 1 + S 2 = 2>, 
2,S 2 = 0,S 1 + S 2 = 2>, 
S 2 = 2,S 1 + S 2 = 2>, 



Ex 


> = 


\Ni 


= 2,iV 2 


= 2, Si 


E 2 


> = 


|JVi 


= 2,7V 2 


= 2, Si 


E3 


> = 


|JVi 


= 2,7V 2 


= 2,Sx 


Ei 


> = 


\Nx 


-2,7V 2 


= 2, Si 


E 5 


> = 


\Ni 


= 2,iV 2 


= 2, Si 


E 6 


> = 


\Nx 


= 2,7V 2 


= 2,5! 


E 7 


> = 


\Nt 


= 2,iV 2 


= 2, Si 


E s 


> = 


|JVi 


= 2,iV 2 


= 2, Sx 



S 2 = 2, S\ 
S 2 = 2, Si 
S 2 = 2, Si 



S 2 = 1 >, 
S 2 = 3 >, 
S 2 = 4 > . 



(29) 



Hopping term in JSJ conserves total spin; therefore, energy in each subspace doesn't depend on the z component of 
Si + S 2 , and states \E e >, \E 7 >, and \E S > form orthogonal subspaces that do not mix with any other states. We 
can then use formula analogous to (|18|) to calculate corrections to the energies of these states in the second order in 
t (see Appendix 151 for details): 



(7 



t 2 

6E/ 2 -4— , 

t 2 

6(7 2 -4— , 



6U 2 



(30) 



In (|3U[) we used U 2 << Uo and neglected U 2 relative to Uq in denominators of the exchange terms. 

Boson tunnelling can connect two subspaces, \Ei > and \E 2 >, and three subspaces, \E$ >, \E^ >, and \E§ > (only 
states with the same component of S z have tunnelling matrix elements). Thus, the energies and eigenstates should 
be found by diagonalizing the matrix 



7~l(3 a — E 0a 5 a p — t 2 ^ 



< f3\(a\ p a 2p + a\ p ai p )\m >< m\(a\ p a 2p + al p ai p )\a > 



E n 



Eoa 



(31) 



Here EQ a is the energy of the state a in the zeroth order in t: E01 = 0, Eq 2 = 6C/ 2 ,i?o3 = 3J7 2 ,£'o4 = 3[/ 2 , -E05 = QU 2 . 
Summation over m covers intermediate states that have 1 and 3 bosons in the two wells, have the same total spin as 
states a and /?, and have spins in individual wells that differ from a and f3 by ±1. The explicit calculation presented 
in Appendix B gives the matrix H.p a and its eigenvalues. Energies of the states with total spin zero are 



h = 3C/ 2 



e 2 = W 2 -- 



(32) 



Energies of the states with S = 2 are 



h = 3C/ 2 -4— , 



e 4 = - 9U 2 - 12t* 



es = 7: 9C/ 2 - 12i 2 



t 4 



i 2 



1-44^+40— U 2 + 9U$ ) , 



lU^+i0^rU 2 + 9Ui 
<->0 Uq 



(33) 
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FIG. 4: Eigenstates of the effective spin Hamiltonian for a two site problem with two atoms per site. Energy and t are 
measured in units of Uo, and we assumed U2 = 0.04C/o- The lowest energy states correspond to total spin 5 = (el), 5 = 2 
(e5), and 5 = 4 (e8). 

and each of these states is fivefold degenerate. Energies ei-es are shown in FIG.Q] 

For U 2 > the lowest energy state is a total spin singlet that has some mixture of S = 2 states in individual wells 
when t is nonzero. The next favorable state has total spin 2, I-E5). When the value of i is increasing, the ferromagnetic 
state \E$) becomes the third low lying state. At this point, we have solved the problem for two wells, taking into 
account competition between hopping and the U 2 term (overall Hilbert space for two wells is 36 dimensional). 

B. Effective Spin Hamiltonian for an Optical Lattice 

In the previous subsection we used perturbation theory in tunnelling t to study the problem of two sites with two 
atoms in each well. If we label the two sites 1 and 2, in the second order in t the effective Hamiltonian can be written 
as 

H 12 = 3U 2 [P(Sx = 2) + P(S 2 = 2)] + |q >! |/3 > 2 J a>0 . n<s < 7 |a < S\ 2 , (34) 

Here P(S{i, 2 } = 2) are projection operators into states with spin S = 2 on sites 1 and 2 and J a ,[3-j,s gives exchange 
interactions that arise from virtual tunnelling processes into states with particle numbers {n\ — \,n 2 = 3) and 
(ni = 3, n 2 = 1). The second term of (|34|) includes all initial states (|7)i and \8) 2 for sites 1 and 2, respectively) and 
all final states (|a)i and \P) 2 )- 

Generalization of the effective spin Hamiltonian l|34|) for the case of optical lattice is obviously 

H = 3U 2 J2 p (S t = 2) + J2 \ a >* \P >j J ^,s < lU < S\j, (35) 

i <ij> 

This Hamiltonian is linear in U 2 and, therefore, can be written as a sum of the bond terms 

H = ^Hij, 

<«) 

Hij = ^[P(S i = 2) + P(S j = 2)} 

+ \a >i |/3 >j Ja,/3; 7 ,5 < j\i < S\j. (36) 

Individual terms Hij differ from (|34ll only by rescaling U 2 — > — , where z is the coordination number of the lattice. 
We did not give explicit expressions for J a ,i3 : -y,s in the basis of eigenstates of individual spins £{ and Sj but in the 
basis of eigenstates of the total spin of the pair l|29|l . expressions for J a ,p;-y,s can be obtained from eigenstates and 
eigenvalues of (|34|l (see equations l|30ll - l|33|) and Appendix B)(with a rescaled U 2 ). Therefore, we can write 

E \E%,S?>Hce<B%,S°l (37) 

a,P,S5 
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where states \E%j , S" > have been defined in equations l|2"§|l . Expressing \E*j> , S z > via states |iV, = 2,Si = {0, 2}, Si Z = 
—Si... Si > using known Clebsch-Gordon coefficients 

\E%,S Z >= C I : 'M ;S?,s° Wi = 2, S?, Sg > \Nj = 2, Sf, S% >, 
we can write the Hamiltonian H35f) as 

H = ^2 \ a >i \l >i W Q) /3; 7 ,<5 < /8|i < S\j, (38) 

where states \a) - \S) belong to the set {S = 0}, {S — 2, S z — —2, 2}, and H a! p ;J j is given by proper rotation of 



C. Phase Diagram from the Mean-Field Calculation 

In this section we study the phase diagram of the system described by the Hamiltonian H38[) using translational 
invariant variational wave functions. Such mean-field approach gives correct ferromagnetic and antifcrromagnetic 
states for Heisenberg hamiltonians in d > 2, so we expect it to be applicable in our case. We think that this approach 
successfully captures the main features of the system: first order transition between the spin gapped and the nematic 
phases, the nature of the order parameter in the nematic phase, and elementary excitations in both phases. However, 
we cannot rule out the possibility of more exotic phases that fall outside of our variational wave functions, for example, 
the dimerized phase discussed in |3l| , and numerical calculations are required to study if such phases will actually be 
present. 

As we saw in the previous chapter, energy of the two-well problem is minimized when total spin is 0. However, 
energy on all bonds cannot be minimized simultaneously, so wc cannot solve a problem exactly for a lattice. We use 
a mean field approach to overcome this difficulty, taking variational wave function 

\y>=Y[(c Ofi \N l = 2,S = 1 S z = Q>+ c 2 , m \Ni = 2,S = 2,S z = m>), (39) 

i to=-2,...,2 



|c ,o| 2 + l c 2,™| 2 = 1 - (40) 

m=-2,...,2 

Now we can evaluate expectation value of energy over variational state l|39|) and find the ground state numerically. 
We parameterize (PU)l as 

co,o = cos 6, c 2 , m = sin 9a m , (41) 



Kl 2 = i. 

m=— 2,...,2 

In Appendix |0 we demonstrate that for a region of 9 where the mean-field energy is minimized, [a m ] has the form, 
up to SU{2) rotations 

[am] = (0,0,1,0,0) T . (42) 
Mean field energy does not depend on d and we find in the region of interest the energy per lattice site to be 

zt 2 

E[0] = 3U 2 sin 2 6 + -——(-51 + 4cos26» + 7cos46> - 8^2 sin 26 + 4%/2sin46»). (43) 
12(7o 

One can immediately see that if we try to expand this expression near 8 = 0, there is no linear term, but second 
and third order terms are present, which indicates that by changing the parameters of the Hamiltonian, we will have 
a first order quantum phase transition, at which the value of 9 that minimizes the energy changes discontinuously. 
This is typical for ordinary nematics|42j since in Landau expansion third order terms arc not forbidden by d — > —d 
symmetry. The reason why our transition is of the first order can be traced back to the fact that mean field energy 
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FIG. 5: Dependence of the energy functional l)4M^ on 8 when zt 2 /(U0U2) ~ 0.4928 (energy is given per lattice site in units of 
U a ). 

has terms which mix co : o and C2, m in odd powers, i.e. (fy 002,-2(02 -i) 2 j an d overall U(l) symmetry doesn't prohibit 
odd powers of 9 in H43[) . 

Since phase transition is of the first order, transition is characterized by several regimes. First, when t is small, 
global energy minimum is at t = and there are no other local minima, i.e. we have spin singlets in all individual 
sites. Then, when condition zt 2 _/ (U0U2) ~ 0.4928 is satisfied, a local minimum appears at 0- f» 0.25, see Fig. [S] 
As we continue increasing t, the minimum at nonzero 9 becomes deeper, and eventually at zt^j {XJ§U2) — 1/2 the 
global minimum of l|43|l is reached for sm9 c = 1/3, (see FigJBJ. However, there is still a local minimum at 9 = 0. If 
we keep increasing t, the minimum at 9 = becomes completely unstable at zA. / (U0U2) — 9/16 and there is only 
one minimum at 9+ ps 0.5 (see Fig. 0). As we increase t further, sin# + continues to grow, approaching the value 
sin = (2/3) 1 / 2 . It is useful to point out that when t is changed in experiments (e.g. by changing the strength of 
the optical potential |43|). we expect that the system will not switch between the singlet and nematic phases at t c , but 
will remain in the appropriate metastable local minimum until it becomes completely unstable. So, in experiments 
with increasing t, the transition from the singlet to the nematic states will occur at t + , and in experiments with 
decreasing t, the transition from the nematic to the singlet state will take place at f_. We note, however, that the 
difference between different t is quite small. 

In Figure [5] we show the phase diagram for the insulating phase with N = 2, including a true first order transition 
line at t c and limits of metastability at t_ and t + . The Supcrfluid-Insulator transition line is shown as an eye guide for 
the case of small enough U2/U0; its exact position will be presented in section IVT1 It is useful to point out that in the 
discussion above we used canonical ensemble (fixed number of particles) rather than grand canonical ensemble (fixed 
chemical potential) to discuss the singlet to nematic phase transition. However, intermediate states that contribute to 
exchange interactions always involve one particle and one hole. Hence, their energy does not depend on the chemical 
potential. This explains why the singlet to nematic phase boundary in Figure [5] does not depend on \i. It is consistent 
with our physical intuition that insulating states have a certain number of particles, but their chemical potential \x is 
not well defined as long as /1 is inside the Mott gap. In the discussion presented in this section, we assumed that the 
system remains deep in the insulating phase and the superfluid to insulator transition does not preempt the isotropic 
to nematic transition inside the insulating lobe. Precise conditions under which this is justified will be given in Section 

ED 



D. Quantum Fluctuations Corrections for the Spin Singlet State 



For small enough t, mean-field analysis of the previous section predicts the singlet ground state that does not 
depend on t. Now we will consider quantum fluctuations near this state to obtain more accurate wave function and 
excitation spectra. We can rewrite (|38H via Hubbard operators 

Af = \a)Mr- (44) 

Here \a)i and \/3)i belong to the set {S = 0}, {S = 2, S z = —2, 2}. Commutation relations between Af 13 are very 
simple: 

[Af,A2 s ]=5 fj ,Af-5 aS Af. (45) 
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FIG. 6: Dependence of the energy functional l|43^ on 6 when zt 2 / (U0U2) = 1/2 (energy is given per lattice site per lattice site 
in units of Ui~). 
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FIG. 7: Dependence of the energy functional 143H on 9 when zt 2 / (U0U2) = 9/16 = 0.5625 (energy is given per lattice site in 
units of L r 2). 

Now we introduce boson operators b ai that create states with Si — 2,Si Z — a, and c\ that creates a singlet on ith 
site. Our physical subspace is smaller than generic Fock space of these bosons and should satisfy the condition 

4*+J2 b l b *°< = 1 ( 46 ) 

a 

on each site. 

One can easily check that if we set A a P = b\fip for spin S — 2 states and similar substitution with c bosons when 
one of the states is a singlet state(which we will denote as s), then commutation relations (|45() are satisfied. Since 
for small enough t only a singlet state is occupied in mean field approximation, we can resolve constraint (1461) using 
analog of Holstein-Primakoff representation near Sc = 1 state which is given by 

Af = b\ a b w ,Ar = (1 - b\ b ifi f/% ia , (47) 
A™ = 1 - b%b Wl Ar = bl(l b\ b l0 fl\ (48) 

Now we expand our initial Hamiltonian in terms of now independent operators b\ up to the second order: 

Jj(2) _ ^ (Hap^sb^bfo + H SS ; a pb a ibf3j + Has-spb^bpj + Hsatfsbljbpi) + 
<ij> 

+ 7} ^(Has-fisblfipi + H sa , s pbl i bf 3i + H ss , ss (l - b ] ai b ai - b ] aj b aj )). (49) 

i 

Calculation of matrices H a p n s gives necessary matrix elements 
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We rewrite our Hamiltonian in terms of Fourier transforms of b ak , b a k- 

H {2) = E + | ^ lk{H a p. ss b^ ak b\_ k + H sa . ia pb a kbf3-k + H as - tS pbl k bpk + H sa -p s b\ k bp k ) + 

k 

+2(H - H ss . ss )b ] ak b ak , (50) 

where Eq is classical energy and 

Now we will use canonical Bogoliubov transformations to diagonalize this Hamiltonian. Since most of the terms are 
diagonal in a, (3 subspace, it is easy to see that required transformation mixes operators {b ok ,bo-k}, &-l-fc}j 

{bl lk ,h-k}, {bl k ,b- 2 -k}, and {&L 2fc , 6 2 _ fc }. 
Transformation that mixes the first pair is 

b ok = cosh6> fc /3ofc + sinh0fc/3j_ fe , 



b -k = cosh6> fe /? -fc + smh0 fc /jj fe , 

and complex conjugates. Substituting this transformation into (|50|l and requiring that terms with P\ k !3\_ k and 
(3okPo-k vanish, we obtain the equation for 9k 

tanh20 fe = -^ = ~ 3 ^ lk 



where 



8zt 2 



9k = oTrTfc. 

6 Uo 



fk = 3J7 2 - oTrT fc - 
6 Uo 

Energy of this excitation becomes 



E(k) = sJJt gl = J U 2 (9U 2 (52) 

Equation l|52l) suggests that the first instability appears at k = and gives the phase boundary that agrees with 
the metastability line i+ found in the previous subsection. However, results of the previous subsection suggest that 
phase transition is of the first order and takes place before the mode softening at k = 0. The first order transition 
may also be obtained with the formalism presented in this section by noting that expansion of (|48|l allows third order 
terms ctftf + c.c. 

We can use similar analysis to discuss excitations with other spin quantum numbers. For example, excitations 
with S z = {+1, —1} are diagonalized by analogous Bogoliubov transformations with 6k — > — 6k, and excitations with 
S z = {+2,-2} are diagonalized with transformations with the same 6k- As required by the spin symmetry of the 
singlet state, all of these excitations have the same energy. 

Now we can discuss the approximations made while expanding over b ak> bpk- While transformation (|48|) is the exact 
resolution of the constraint (|46|) . expansion to the second order adds states with higher boson occupation numbers 
and changes Hilbert space (this is completely analogous to usual antiferromagnet spin- wave theory). However, if a 
posteriori we can verify that only states with occupation numbers n a i = {0, 1} are present in the ground state, then 
expansion of the constraint H46fl up to the second order was justified. The parameter that controls such expansion is 

b^bai = ^J2 b ik b a k = y"sinh 2 9kj^s- 

Calculation of this quantity while the singlet state is still a global maximum for d = 3 gives numerical values < 0.001; 
therefore, our expansion is much more precise than for Heisenberg antiferromagnet, where this quantity is not much 
smaller than 1 and one needs the condition S 3> 1 to justify the spin wave theory. 
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FIG. 10: Dependance of excitation gap (measured in units of U2) on t, measured in units of 
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E. Spin Wave Excitations in the Nematic Phase 

Now we will consider excitations for the states with nematic order. For the states described by l|39|) - (|42J) . there 
are no expectation values of the spin operators (S) — 0, but there is a nematic order l|27l) when 9 7^ 0. For example, 
when d is pointing along z we find 



Q ab = sin 2 9 



( 1 
1 

v -2 



(53) 



In the singlet phase, expectation values of both the spin operators are (S) = and (Q a b) = 0. In the singlet phase the 
system has a gap to all excitations of order U2, while nematic phases have gapless spin- wave excitations that originate 
from the continuous symmetry breaking. The general form of the state with minimum energy is expressed via Euler 
angles of order parameter d as 

^( 7 )C/ y (/?)C4(a){0,0,l,0,0} T , 

where ^(7) are finite angle rotation matrices. From 15311 and i|28[) we can express the nematic order parameter for 
such a state as 

Qab = -3 sin 2 9(d a d b - -S ab ). 

Goldstone theorem tells us that low lying modes will be fluctuations of direction of d, and there will be two degenerate 
modes. We can utilize the approach used in the previous subsection to consider excitations in the nematic phase. 
Here, we should make generalized Holstein-Primakoff expansion near the nematic state. First, we make unitary 
transformation in Hilbert subspace of each site, which is given by 



V 





/ 


cos 6 


1) 







2) 







3) 




— sin 


4) 







5}/ 


V 






sin(9 \ 

1/V2 1/V2 
-1/V2 1/V2 

cos6» 

1/V2 000 1/V2 
-1/V2 000 1/V2J 



( \S = 0,5 = 0) \ 
S = 2,S Z = -2) 
S = 2,S t = -l) 
\S = 2,S Z = 0) 
|S = 2,S, = 1) 
V \S = 2,S Z = 2) J 



(54) 



Making appropriate transformation on H a p. 1 g 1 we can write our Hamiltonian as 

<»i> 



(55) 



where states \a) - \S) belong to the set {|0) — |5)}. After that, we proceed exactly as in the previous subsection, 
expanding near |0) state. Since dependance 9 on t 2 /(U0U2) is determined by the minimization of the energy, linear 
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FIG. 12: Energy of the gapped excitations in the nematic phase of N — 2 atoms per site at zero wavevector. The energy gaps 
and t are measured in units of U2 and \ I U ° U2 , respectively. 



terms in bka and b k are absent. Quadratic terms have exactly the same form as in 149f) . and all matrices become 
diagonal due to the proper basis choice l|54|) . Now we can use Bogoliubov transformation to diagonalize the quadratic 
part. For excitations to states |1) and |2), we obtain energy dependance 

Ef(k) = E 2 2 (k) = i(-16 7 ^ (4cos20 + x/2sin20) 2 + 

(877- - 18 7 fe^ + 9C/ 2 + (2( 7fc - l)^- + 9U 2 ) cos 26 - 7^ cos 40 
Uo Uq Uq Uq 



+4V2(1 - lk )%- sin 26 - 4^2-^ sin40) 2 ), 

where 7^. was defined in i|51[l . and dependance of 6 on zt 2 / (U0II2) is shown in Fig. [SJ We find that for k — 0, energies 
of these excitations are zero, as expected for nematic waves from Goldstone theorem. These excitations create states 
with S z = ±1. For small k, the energy of excitations depends linearly on |fc|, and dependance of spin wave velocity 
on the parameters of the lattice is shown in Fig. ^2 

Let us now consider gapped excitations for the nematic phase. Excitation to the state |3) corresponds to longitudinal 
fluctuations in the value of 6, and the energy of such excitations becomes zero at t- since at this point fluctuations 
of 6 are not suppressed. Excitations to the states |4) and |5) correspond to the creation of S z = ±2 states and they 
are degenerate. For all of these excitations, energies are minimized for k = 0. Dependence of the gap on parameters 
is shown in Fig. 
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F. Effects of Small Magnetic Field 



Let us now consider the effect of a small magnetic field, [iH << U%, on our system. For U2 in the range of kHz 
(see Fig. |5J), this corresponds to magnetic fields smaller than a lmG. We suppose that the field is small enough that 
it does not change the scattering lengths due to the energy level shifts inside of atoms. Since all atoms have the same 
gyromagnetic ratio, interaction with external magnetic field depends only on the total spin, and the internal structure 
of the states is not important. In the case of a nematically ordered insulating state, the ground state energy doesn't 
have any contributions linear H (this follows from < 0|HS|0 >= 0), and the second order contribution depends on 
the relative orientation of the nematic order parameter d and magnetic field H . Suppose d is directing along the z 
axis and H lies in x, z plane. In the second order perturbation theory, energy correction to the ground state is always 
non positive: 

zp(21 < 0|/iHS|n >< n|^HS|0 > v^. 2 < 0\S x \n >< n\S x \0 > 

tr = — — 

this quantity is of order — (/xif x ) 2 /t 2 . Since for the state with d\\H it is again zero and the system doesn't benefit 
from magnetic field, energy is minimized when d lies in a plane perpendicular to i?(this is completely analogous 
to Antiferromagnet). Using this property, one can a distinguish nematic phase from a singlet phase. One should 
apply a small magnetic field in z direction to fix the plane in which d lies, release the trap, and let the atoms fall 
in the gravitational field with some magnetic gradient, to separate the states with different S z . Then, one should 
measure quantities of each spin component. These values will have a sharp change when we cross the first order phase 
transition line. Knowing how to express spin states via original boson operators, we can calculate expectation values 
of different spin components to be: 

m = n-i = I cos[6>] 2 + cos[0] sin[0] + j sin[0] 2 , 

00 D 




FIG. 13: Dependence of occupation numbers no and m = n_i on t for an insulating state with two bosons per site, N = 2. 
Tunnelling t is measured in units of ■ ' u ° 1,2 
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V. LARGE NUMBER OF PARTICLES PER SITE 



In this section we discuss the case N ^> 1 for both parities of N. We show how one can separate variables describing 
angular momentum and the number of particles in each well|44|. and derive an effective Hamiltonian which is valid 
under conditions U2, Nt <C Uq, which is less restrictive than in section Hill for N >> 1. 

When we have N spin-one bosons localized in a well in the same orbital state, their total spin may take any value 
that satisfies constraints 

S + N = even, (56) 
S < N. (57) 

We define 'pure condensate' wave functions as 

\N,n) = ^(n x 4 + % 4+ n z 4) N \Q) (58) 

that minimize the U2 interaction energy at the Gross-Pitaevskii (mean-field) level in each given well [2]]]. Here 
Af = [2(N — l)!] 1 / 2 is a normalization factor, which is calculated in Appendix iDl 
Now we can construct states as 



N 



V(n)|iV,n), (59) 



where J n stands for J da/ (Air). 

Condition (|56|) corresponds to the symmetry of the states (|58|l 



\N,-fi) = (-) N \N,n). (60) 



Hence, we need to consider only wave functions that satisfy ip(—n) = (—) N ip(n). 
Now we can consider how a spin rotation operator acts on the wave function ^>(n): 



)N 



^(n)e- ie ° s °|iV,n) = / i>{e w " s "n!)\N , n'> = 

in' 

■tp(n~t + e a p y 6 a n'p))\N,n'). (61) 



Expanding the last expression for small 8 we find 

d 

L a tp = -teapynp-g—ip, (62) 

where we used L a rather than S a to show that it acts on the wave function ip. Therefore, operator L is an angular 
momentum operator for n. If we want to construct any spin state, we should take tp(n) to be usual spherical harmonic. 
We note that the S = result in |2l| is just a special case of our general statement. The most general form of the 
state in a well can be expanded as 

= ^2 c L>m Y Lm (n), 

\m\<L, 

where L satisfies conditions l(55 )l -(|5T )l . 

At this point, what we have done is valid for all N, not only big ones. This representation is particularly suitable 
for N ^> 1 since in this limit states that correspond to different n's are orthogonal to each other(see Appendix IDl) 

(Ar,n 1 |Ar,n 2 ) = ^(n 1 -n 2 ) (63) 

The delta function is defined from the condition 

I [ fh(n 1 )ff f (n 2 )8 N (n 1 -n 2 )= [ fh(n)f N (n) (64) 

• ' ri j •> 112 •' 11 

for the functions that satisfy /n(— n ) = ( — ) iV /iv(n). 
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We show in the Appendix iDl that 

al\N,n) = (JV + ljVa^jJV + l.n), 

a a \N,n) = N 1/2 n a \N - l,n) (65) 



after projecting into the 'pure condensate' wave functions. 

This allows us to represent (|65|1 as the product of two operators, which act in different spaces. For each trap we 
define the particle creation and annihilation operators that change the number of particles N but not the direction 
of n[H 

bllN,,^) = (^ + l) 1 / 2 |iV l + l,n J ), 

b t \N t ,n % ) = Nt /2 \Ni -l,ni). (66) 
The number of particles in each trap may be expressed using b operators as 

Ni=b\bi. (67) 

Hamiltonian ijjjj can now be represented as 

H = -t n 4 n, (bib, + b]bi) -^^i + Y^^^-^ + Y - 2 ^)' ( 68 ) 

(ij),a i i i 

where 

d 

Li =-m,x - — . (69) 
drii 

Now, if we are in the Mott insulating phase, we can easily derive the effective Hamiltonian for ?/>(n). Using the 
second order perturbation theory, we hnd the effective Hamiltonian on the sphere to be 

U 2 ^ f2 2NH 2 x . 
£ Uq '-rr 

V(n) = (-lfi;(~n). (70) 
We note that this Hamiltonian corresponds to a lattice of quantum rotors that interact via quadrupolar moments. 



A. Mean Field Solution 



Now we can find a mean field ground state of H70|) . We consider the case when the quadrupolar interaction term 
is much bigger than the kinetic term. We show that in this case the ground state is a uniaxial nematic and find its 
energy. Comparing the energy of this state to that of a singlet, we estimate the phase boundary for nematic-singlet 
transition for even N. 

Our general mean field anzats has the form 

i*>=n*(fc.p*)> / i*(mi 2 =i- ( 7i ) 

Expectation value of i|70|) per well over the wave function l|7f|) equals 

— < L i > — J < n a np >< n n a >, 

where 
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-9 1 d , d , 1 d 2 

sin 9 89 86 sin o^ 2 

As in the case of N = 2, all of the states that can be transformed into each other by global rotation have the same 
energy. Therefore, we can impose three additional conditions. The best choice is to require symmetric real matrix 
< UiTij > to be diagonal and to choose < n 2 > to be the biggest eigenvalue. In such a gauge interaction, the term 
becomes — J(< n 2 > 2 + < n 2 > 2 + < n 2 > 2 ). Since we have the extra constraint < n 2 > + < n 2 > + < n 2 z >— 1, it 
is now obvious that interaction energy is minimized when 

< n 2 >= cos 2 6 — * 1. 

However, states with sin 2 # — ► have higher angular moments, and the ground state is determined by the competition 
of these two factors. We write mean field Gross-Pitaevskii equations to determine the ground state 

j U 2 1 d . d 1 d 2 



2 sin 6*961' 09 ' sin 2 9 dip 2 



2J(n 2 < nl 



> +n 2 y <n 2 y > 



-nl<nl >)}*(0,<p) = X9(0,<p), 



(74) 



where A is Lagrange multiplier. Now we consider the case J 3> U 2 . In this case interaction energy dominates and we 
expect sin 2 9 — > and, therefore, wave function becomes localized near z and —z directions. We can solve the problem 
by expanding only near 9 = and then taking an (anti)symmetric combination to satisfy {ZD} - It is obvious that if 
we expand the kinetic part of l|74|) up to the first nonvanishing order in 6, then we will get a two dimensional Laplace 
operator A ra:cin , and our problem becomes equivalent to a harmonic oscillator. Effective parameters are expressed as 



1 



= AJU 2 <n 2 z ~ nl y 



> 



(75) 



2 > while obtaining a harmonic hamiltonian from i|74f) . 
~5l, i.e. the ground state is a uniaxial nematic. Since 



Since we have already neglected higher order terms in < 
with the same accuracy we can set u> x — u> y — \fAJlJ2 in 
we know wave functions, we can calculate the expectation value of energy. We will use the fact that for a harmonic 
oscillator the expectation value of kinetic energy is the same as that of potential energy. Energy of the ground state 
becomes 



E 



,)- J< n\ >. 



Quantum fluctuations of the direction of n equals 



< nt >=< nt >= 



(76) 



and expectation value of the ground state energy is 

„ 3 JU 2 
E = -u H 

4 U! 



J = 2y/U2J - J. 



(77) 



Symmetrization or antisymmetrization of the wave function introduces exponentially small shifts in energy, so in limit 
J»£/2 energy doesn't feel the parity of N. Though in this subsection we explicitly started from variational anzats 
l)7ip. now we can justify it in the limit J ^> U 2 since in this case quantum fluctuations of the direction of n are small 
and given by (|76|l . 

From section Hill we know that for small nonzero t there is a uniaxial nematic state for odd N. Since in the opposite 
limit there is also a nematic state, we expect that for all N ^S> 1, N = 2n + 1 the insulating state will be nematic. For 
the case of even N, there is always a singlet state in which mean-field energy equals —J/3. Comparing this energy 
with H77|) . we can estimate the first order transition point as U 2 = J/9. At this point 



< n\ >= 1/12, 



so we expect our expansion to be valid. 
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VI. GLOBAL PHASE DIAGRAM 

In the earlier sections we have established spin structure of insulating phases of S — 1 bosons in the optical lattice in 
various limits. Here we summarize our arguments and discuss implications of our results for the global phase diagram. 

A. Two and three dimensional lattices 

In two and three dimensional lattices, insulating states with one atom per site are nematic as long as the perturbation 
theory approach in t/Uo remains valid, as was shown in section ITTT1 For an arbitrary odd number of particles per site, 
N, and in the limit of small tunnelling (Nt) 2 /U~o << U 2 , the nematic order in the ground state was also established 
in section [ITTI For large, odd N, the nematic order in the ground state can be proven when (Nt) 2 /Uo becomes larger 
than U2 (but still smaller than Uq), as was demonstrated in section It is also natural to expect that the superfluid 
polar phase develops from the nematic insulator (both states break spin rotational symmetry without breaking the 
time reversal symmetry), so we expect the nematic order even when Nt/Uo is not small and the system is close to 
the superfluid-insulator transition. In all cases we find that insulating phases with an odd number of particles per 
site are nematic. 

In the case of two particles per site, the results of section llVl establish that for small enough U2/U0 there is a first 
order transition between the spin singlet phase (for small t) and the spin nematic phase (for larger t) at zt 2 /UoU2 = 0.5 
(z is the coordination number of the lattice). Analogously, for large, even N, results of sectionlVlshow that the singlet 
insulating ground state goes into spin nematic at zN 2 t 2 /U~oU2 = 9. Since for small enough U2/U0 we expect nematic 
spin order close to the SI transition into the polar superfluid phase, we propose that in this case insulating phases 
with an even number of particles per site are either singlet or nematic with the first order transition at some critical 
value of tunnelling t c . 

In all of our earlier discussions, we assumed that Mott insulating lobes for even fillings are big enough to have the 
transition into the nematic phase before superfluidity sets in. This assumption is controlled by the smallness of the 
ratio U2/U0. Here we will discuss the superfluid - insulator phase boundaries and estimate how small U2/U0 should 
be for the singlet-nematic transition to lie inside the Mott phase. 

Assuming transition from the spin-singlet insulating phase, the mean-field calculation of the superfluid-insulator 
phase boundary was given in |26j | . Analysis presented in this paper shows that the critical value of tunnelling, after 
which the Mott phase doesn't exist, is given by 



Uq + 2U 2 



-{2N + 3 + 2^N 2 +3N). 



zt S i 3 

We will use this critical value tsi as an estimate of the superfluid-insulator transition. For N — 2, singlet-nematic 
phase transition takes place at zt 2 /UoU 2 = 0.5. The condition t c < tsi for N — 2 is satisfied, if 

For the case N >> 1 the requirement of t c < tsi becomes even more restrictive, namely 

zU 2 

One can see that depending on the exact value of ZU2/U0, there are different possibilities for Mott lobes with an even 
number of particles. When ZU2/U0 < 0.01, all insulating phases with even filling factors are spin singlet for small 
tunnelling and spin nematic for larger tunnelling. For 0.01 < ZU2/U0 < 0.1, insulating phases with small, even filling 
factors have both singlet and nematic regimes, but insulating states with sufficiently large even fillings have only the 
singlet phase. Finally, for 0.1 < ZU2/UQ, all insulating phases with even filling factors are in the spin singlet state. In 
Figs. 1141 and [TSI we combine these results with the schematic representation of the SI transitions to obtain the global 
phase diagram. 

B. One dimensional lattices 



For one dimensional lattices we established that when N 2 t 2 /Uo << U2 the system will be in a uniform singlet phase 
for even fillings and in a dimerized singlet phase for odd fillings (when there is only one atom per site the dimerized 
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FIG. 14: Global phase diagram for S = 1 bosons in 2d and 3d optical lattice, for ZU2/U0 < 0.1. Mean-field analysis of 
the superfluid to insulator transition was done in |2?|. In this work we concentrated on discussing the spin structure of the 
insulating lobes. Singlet-nematic first order phase transition for N = 2 takes place for zt 2 /U0U2 = 0.5 (z is the coordination 
number of the lattice). For large, even N singlet-nematic phase transition occurs at zN 2 t 2 /U0U2 = 9. t c marks the actual first 
order phase transition and t- and t+ are the limits of metastability. Note that the system may also have fragmented superfluid 
phases for small t that are not shown here 
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FIG. 15: Global phase diagram for S = 1 bosons in 2d and 3d optical lattice, for zUi/Uo > 0.1. Superfluid-Insulator transition 
for even filling factors takes place before singlet-nematic transition. 



phase has been verified in the regime t << Uq). The nature of magnetic order close to the tips of the insulating lobes 
(when the perturbation theory in t is not applicable) is less clear. However, we expect that the phase diagram for 
the one dimensional lattice is qualitatively similar to two and three dimensional cases with one important difference: 
instead of the nematic phase, we have dimerized singlet states. This will be discussed in future publications. 
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FIG. 16: Probing dispersion relation using Bragg scattering. 



VII. DETECTION OF SPIN ORDER IN INSULATING PHASES 

Now we discuss two approaches to detection of the singlet-nematic phase transition for S = 1 bosons in an optical 



lattice. One way of detecting such a transition has already been noted in section HV Fl where we proposed to introduce 
an easy plane for nematic order by applying a small magnetic field, then releasing the trap and measuring the number 
of particles of different spin components. Spatial separation of different spin components can be achieved by applying 
magnetic field gradients during the free fall of the atoms. For the case of N = 2, with a small magnetic field applied in 
the z direction, expectation values of n(S z = 0) and n(S z = 1) = n(S z = —1) have been calculated and are shown in 
Fig. Since the phase transition is of the first order, there is a sharp change which can be measured experimentally. 
We note that N = 2 case also have particular experimental advantage over other filling factors due to the absence of 
three particle losses and the least restrictive condition on U2/U0 for observation of singlet-nematic transition. 

The second approach to experimental detection of singlet and nematic insulating phases relies on the measurement 
of excitation spectra. As discussed in sections IIV Dl and IIV El the singlet phase has a nonzero gap to all excita- 
tions, whereas the nematic phase has gapless spin wave excitations. To measure the excitation spectra, we propose 
using Bragg spectroscopy, which was used successfully to identify sound-like Bogoliubov excitations in condensates 
of spinless particles |45j. In such experiments the optical lattice should be illuminated by two laser beams with 
wave vectors ki and k2 and a frequency difference uj, which is much smaller than their detuning from an atomic 
resonance. The intersecting beams create a periodic, travelling intensity modulation that creates external potential 
due to ac Stark effect of the form V a /3 cos (qr — tot). Here we introduce spin indices for the external potential since 
the ac Stark effect may introduce mixing between different S z components. A response of the system to such po- 
tential may be calculated using Fermi's golden rule. Interaction is expressed in second-quantized notations such as 
Vai3/2(p'l f3 ((i)e~ lujt + p^ /3 (-q)e l " t ), where pl p = J2k a lk+ q a 0k- The scattering rate is given by 



Y E \(f\V a ^\g)\ 2 5(nw -E { + E g) . 
/ 

If the resonance state is far detuned from the excited states, then V a p has the form VS a /3, and couples only to the 
total number of particles in each well and doesn't feel internal spin structure. Low lying excitations in insulating 
phases don't change the number of particles on individual sites, so VSap interaction won't produce any Bragg peaks 
for low lying excitations. Therefore, it is necessary for detection that V a g deviate from VS a p, which can be achieved 
by making detuning comparable to level spacing of fine and hyperfine components. From section TlV El we know that 
for N = 2, nematic spin wave excitations correspond to S z = ±1, longitudinal excitation corresponds to S z = Q, and 
there are also gapped excitations with S z = ±2. Since the nematic state has S z = 0, it is necessary to have nonzero 
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FIG. 17: Positions of Bragg peaks for |q| = 0.02, z = 6 in the nematic phase with N = 2 atoms per site. Energy is measured 



^b,±i 5 ^b.cb and V±i i=F i to observe each kind of these excitations. In Fig. El we show dependence of N — 2 peak 
positions on t for fixed q. 

Finally, we consider the effect of inhomogeneous trapping potential. When this local trapping potential tpi varies 
smoothly from site to site, it is not the chemical potential //$ which is fixed across the trap, but the sum ip^ + /!;. 
Therefore, if ifi varies considerably, we will have insulating regions with different occupation numbers as well as 
regions with superfluid order, all in the same trap, as was discussed in for the case of spinless atoms. Therefore, 
Bragg scattering experiments for fixed q will exhibit resonances coming from the regions of the lattice at different 
filling factors. Relative intensity of these resonances will be determined by the relative number of particles in each 
region. Interpretation of Bragg experiments will be easier if the trapping potential is not harmonic, but has sharp 
borders, so the whole system has essentially the same density. 



In summary, we have considered Mott insulating phases of spin-1 atoms with antiferromagnetic interactions in 
optical lattices. In the experimentally interesting limit U2 <C Uq, an d deep inside the Mott phases Nt <C (N is 
the filling factor), we performed detailed calculations for the following cases: i) odd number of particles per site and 
(Nt) 2 /Uq <C U2; ii) two particles per site and an arbitrary ratio of t 2 /Uq and U2; iii) large number of particles per 
site N ^> 1 with an arbitrary ratio of (Nt) 2 /Uq and 1/2- Based on this analysis we argued that in two and three 
dimensional lattices insulating phases with an odd number of particles per site are always nematic. For an even 
number of particles per site, there is either a spin singlet phase or a first order phase transition between spin singlet 
and nematic phases controlled by the depth of the optical lattice. The resulting global phase diagrams are shown in 
Fig. ED an dEl We have considered excitations for singlet and nematic phases and have reviewed the effects of small 
magnetic field. For one dimensional lattices we have found dimerized singlet phases for insulating states with odd 
fillings. We also discussed different experimental techniques to identify the proposed phases. 

We thank E. Altman, D. Haldane, W. Hofstetter, D. Podolsky, S. Sachdev, A. Sorensen, DW. Wang, and F. Zhou 
for useful discussions. This work was partially supported by the NSF Career Award DMR-0132874, PHY-0134776, 
and by the Sloan and Packard Foundations. When this work neared completion, we learned that similar results have 
been obtained by M. Snoek and F. Zhou|4r|. 

APPENDIX A: DERIVATION OF THE EFFECTIVE MAGNETIC HAMILTONIAN FOR AN 
INSULATING STATE WITH ODD NUMBER OF ATOMS 

To be able to derive Jo, J\, J2 dependance on n, we should know how to write down explicitly, in terms of creation 
and annihilation operators, all of the states that we are interested in. To do this we introduce singlet pair creation 
operator (summation over repeated indices is presumed over {x,y,z}): 
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(alf - 2ala\ 
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which has the following commutation relations: 

[a p , 6*] = 2al [a , fit] = 2a\, [a ± ,9^] = -2a%, 

[a^tf] = 0,[S a ,tf] = 0,[S 2 ,tf] = 0. (Al) 

Since 8^ commutes with spin operators, it doesn't change total spin and spin components but just adds 2 electrons. 
Therefore, we can construct unnormalized spin states for arbitrary n in the following way: first, write down a state 
with necessary spin for a small number of particles; second, apply 6^ as many times as needed to get the desired 
number of particles. 

Using this procedure we obtain that states 

aUe t ) n |0>,aL(^ t r|0> ! 4(^) n |0> 

belong to S — 1 low energy subspace and are orthogonal since they are different eigenvectors of S z . Making orthogonal 
transformation that leads to {aj,, a' y , a|} basis, we can write three orthonormal states with S = 1 as 

\S = l,p, 2n+ 1 >= J_ o t(gt)n| >, 

where normalization factor 

(2n + 2s + l)!! 



f(n;s) = s\n\2 T ' 



(2s + 1)!! 



was calculated in |22|. In our calculation later we will need more general normalization factors, so first we will derive 
the way to normalize our spin states. 

1. Normalization of the states 

We will be interested in normalization of the states 

\a,b,n>=alal(9^ n \0>, 

and calculation of 

f(a,b,p,q,n) =< a,b,n\p,q,n >=< 0\{6) n a a a b ala\{0^) n \0 >, 
where a, q S {x,y, z}. Let's consider a coherent state 

e alx 1 +alx 2 +alx 3 ^Q > ^ 

We observe that this state is a linear combination of Fock states, with the coefficients being polynomials in {x\, X2, X3}. 
To extract the weight of the state \a,b,n >, we need to calculate the quantity 

\a,b,n>= [T^ b (x)e a ^ +a t^ +a ^\0 >] x=0 , 

where 

T^ b (x)^V Xa V Xb (A x ) n . 
We use the normalization condition for coherent states(see e.g. |4g) 



to calculate 

" r " ' ' Vl P ;,We"*" '—"J|x=0,y=0 



[TZ b (x)T"{y)e XlVl+X2y2+X3V ^ 
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T n ( , T n ( s (gigi + X 2V2 + ^ 3 ) 2 " +2 
(2n + 2)! 

We can expand T™ b {x)T™ q {y) using the extended Newton binomial formula : 



ni+n2+ri3— n 



V72mi Y72m 2 y72m 3 y7 y 



mi!m 2 !m 3 ! 

mi+m2+rri3— n 



Let's first consider the case when sets {a, b} and {p, q} coincide. We have two essentially different cases: a = b and 
a 7^ b. Without loss of generality, suppose for first case a = b = p = q = x. Then, 



n> 2 1 

— i — i — 



ni!n 2 !n 3 ! hH 2 U 3 \ 

ni+ri2+n3=n,'l+'2+'3=2n+2 



(V :cl V !yi ) 2 ^+ 2 (V :C2 V, 2 ) 2 " 2 (V :C3 V y3 ) 2 ^( a;i2 ; 1 )^( a ; 2 y 2 )^(x32/3) i3 



E ( , ) 2 (2n 1 + 2)!(2n 2 )!(2n 3 )!. 

^ ni!n 2 !n 3 ! 

This double sum can be calculated; the answer is 

2 

f(x, x, x, x, n) = — (3 + 2n)(5 + 3n)(2n + 1)!. 
f 5 

For the case of different indices, the normalization is 

n\ ,o 1 



x 



f(x,y,x,y,n)= V ( — — — -) 



2 X 



, ni!ra 2 !rc 3 ! Zi!Z 2 !Z 3 ! 

ni+n2+n3=n,'l+'2+'3=2n+2 



E ( , ) 2 ^i + l)!(2n 2 + l)!(2n 3 )! = 

ni!n 2 !n 3 ! 

rii+n2+"3=n 

4(3 + 2n)(5 + 2n)(2n+l)!. 
15 

Now let's consider the case when {a, b} and {p, q} don't coincide. There are 4 essentially different cases: 

(x, y, x, z), (x, x, y, z), (x, x, x, y), (x, x, y, y). 

All other normalizations can be obtained from these by proper permutation of indices. In the first three cases, overlap 
is zero since(say, for the first case) a nonzero value comes from the term that satisfies the conditions 

2m + 1 = 2mi + 1, 2n 2 + 1 = 2m 2 , 2n 3 = 2m 3 + I, 

which doesn't have integer solutions. However, in the fourth case, overlap of the states is not zero: it comes from the 
terms obeying 

ni + 1 = mi,n 2 = ra 2 + l,n 3 = m 3 . 
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We can calculate this quantity analogously to the previous calculation, but we can use another trick: 
< O|(0) n+1 (0t) n+1 |O >= 3 < 0\(8) n a x a x at4(8^) n \0 > +6 < O|(0)"a x a x a|,a|,(0t)n|o > . 

Therefore, 

f(x,x,y,y,n) = h(2n + 3)! - 3^(3 + 2n)(5 + 3n)(2n + 1)!) = ^-n(3 + 2n)(2n + 1)!. 
15 15 

Since sometimes it will be more convenient for us to work in {+, — , 0} basis, let's also write down all nonzero overlaps 
in this basis (up to trivial permutations): 

/(±, ±, ±, ±, n) =< O\(6) n a±a±aiai(0t) n \O >= A(5 + 2n)(3 + 2n)(2n + 1)!, 

/(0, 0, 0, 0, n) =< 0\{e) n a o a o alal(9^ n \0 >= A(5 + 3n)(3 + 2n)(2n + 1)!, 

/(0, ±, 0, ±, n) =< 0\(6) n a±aoalal(tf) n \0 >=^-{5 + 2n)(3 + 2n)(2n +1)!, 

15 

/(+, -, +, -,n) =< 0\(9) n a+a-alai(8 r ) n \0 >= 4(5 + 4n)(3 + 2n)(2n + 1)!, 

15 

/(+,-, 0,0, n) =< 0|(fl) n o+o_oJoS(fl t ) n |0 >= -^-n(3 + 2n)(2n+l)\. 

1 5 

2. Calculation of to 

So, first let's calculate the energy for total spin state. From known Clebsch-Gordon coefficients, this state is 

1 



\S = 0, S z = >= — (|1, 1 >i |1, -1 >j -\l,0>i |1, >j +|1, -1 >i |1, 1 >,). 



Rewriting it via creation operators, we get the normalized state 

ISi + Sj = o>= — -at o] ]jel) n (e]) n \o > . 

\/3/(n,l) ' p J 

In the second order of perturbation theory, the energy expectation value is 

t 2 

J2 ' F _ f < m \( a l P a 3P + a] p a ip )\Si + Sj > | 2 . 

The intermediate state |m > cannot correspond to Sj = 2, Sj = since in this case two spins can't add to form a 
singlet. There are always at least two possible states: 

Si = 0, rij = 2n + 2, Sj = 0, rij = 2n, Si + Sj = 

and i <-> j. The matrix element for each of these states is 

-=-^ < O|0 n+1 atat(0t)"|o >,< 0|e n a p aJ(flt)"|o >,■ 1 = 

V3/(n,l) P ?V ' ' ' V(/(«+l,0)/(n,0)) 



' <Or+W(0t r | O >^ M /(n,l) 1 l(f(n+m 



V3f{n,l) p ,v y . . . ^ (/(n+1;0)/K0)) y V 3/(n,0) 
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Finally, this term gives 

2t 2 /(n+l,0)_ 4i 2 (n+ l)(2n + 3) 
~3(£/ -2[/ 2 ) /KO) 3(Uo - 2C/ 2 ) ' 

In general case, n^O, there can also be intermediate states 

Si = 2, n t = 2n + 2, Sj = 2, n 3 = 2n, Sj + Sj = 

and i <-» j because 5, = 2 and Sj = 2 can form a singlet. However, in the case of n = 0, such a term doesn't exist but 
can appear if a,j p , acting on aj g (#j)™|0 >, breaks a singlet part and the result has a mixture of S = 2. But for n = 
there is no singlet part, so these intermediate states are absent. One should also note that even for big n we won't 
have terms with higher spins, i.e. Si = 4, because our perturbation in the Hilbcrt space of i — th site is a vector and 
can have matrix elements only for the states with spins differing by ±1. 

From general considerations we know that for Sj + S } ; = energy has the form 

2 .2/ AM , AM , 
[ U a -2U 2 U + 4U 2 h 

where the term with Uq — 2U2 in the denominator comes from the processes of the first kind, and the term with U0+AU2 
in the denominator comes from the processes of the second kind. We calculated f\ (n) earlier to be 2(n + l)(2n + 3) /3. 
Now, to find f2(n), we take the limit U2 — > 0. In this case we don't need to know the exact form of states with 
Si = Sj = and Si = Sj = 2, Sj + Sj = since their energy is the same. We can take the intermediate state |m > to 
be 

im >= 4^(4, P ^, P + ^i, P )4AM) n ( 6 ]) n \° >- 



where M is normalization of the intermediate state. Then, the second order energy is 

I < m\m > ,A 2 f 2 M 

t/o " C/o3/(n,l)2- (A ' j 

Using formulas from the previous subsection, we can calculate M. Then, we can write an equation 

2(/i(n) + / 2 (n)) = 3/( ^ 1)2 = 2 + ^n(5 + 2n), 

and finally obtain 

4t 2 (n+ l)(2n + 3) 16t 2 n(5 + 2n) 



eo = -- 



3([/ - 2E/2) 15(J7 + 4/7 2 ) 



3. Calculation of ei 



Let's consider the case Sj + Sj = 1. There are no intermediate states Sj = 0, Sj = 2 or Sj = 0, Sj = since these 
states can only form Sj + Sj = 2 and Si + Sj = 0, respectively. Sj = 2, Sj — 2 can add up to form Sj + Sj = 1 , and this 
is the only contribution. Sj + Sj = 1 subspace is three dimensional, and from rotational invariance we can choose any 
state we want to calculate the energy. Let's choose our initial state to be Sj = 1, Sj = 1, Sj + Sj = 1, Sj Z + Sj Z = 0. 
From known Clebsch-Gordon coefficients, we can write this state using creation and annihilation operators as : 



a\ , a\ — a] a\ , , , 

\Si +sj = i >= ,+ ^ /(n X) 3 ' + vtnfy n \o > ■ 

Normalization of 

(4,+ a 3,+ + a l,o a 3,o + 4,- a j,-)(4,+ a i- - 4,-4>,+)\ S i + Sj = l> 
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equals 

2(/(+, 0, n) (-2n) 2 /(+, 0, n - 1) + /(+, +, n)(2n) 2 /(+, +, n - 1)) 
^n(l + 2n)(3 + 2n) 2 (5 + 2n)(2n + l)! 2 . 

Hence, the energy for this case is 

4t 2 n(5 + 2n) 



c-i = - 



5(u + m) ' 



4. Calculation of £2 

Lets choose the normalized state for which we will calculate the second order energy to be the state with total 



Siz + Sj z — 2 



\Si +s 3 = 2 >= i a i + ai + (eine}r\o > . 



/(n,l) 

The intermediate state should also have Si + Sj = 2, Si Z + Sj Z = 2. Such a state can't belong to Si — 0, Sj = 
subspace since this pair of spins can't add up to form Si + Sj = 2. There are four possible intermediate states with 

Si = 2,ni = 2n + 2, Sj = 0, nj = 2n, 

Si = 0,n t = 2n + 2, Sj = 2, nj = 2n, 
and i <-> j. In the first case, matrix element equals 

— i— < 0\e n a + a + aW + (9^) n \0 > t < 0\6 n a + aU6^) n \0 >-• 1 = 

/(n,l) ' V/(".2)/(",0) 



;f{n,2)f(n,l) 



/(n.i) V(/MM VV/(".°)/" 

Finally, this term contributes to the energy 

2t 2 /(n,2) _ 4t 2 (2n + 3)(2n + 5) 
~(^o + ^2)/(n,0) ~ 15(C/ + t/ 2 ) 

The second process contributes with the same energy denominator but with a different dependance on n: 

I6t 2 n(n + 1) 
~ 15(U Q + U 2 ) ' 

Now, let's consider the contribution from the Si — 2, Sj = 2, Si + Sj = 2 intermediate state. Using the same trick 
as at the end of the calculation of eo, we only need to calculate norm of 

+ 4o«i,o + al_aj,-)al + a] !+ (elr(6]r\0 > . 

This quantity equals 

/(+, +, n)(/(0, 0, n - 1) - 2(2n + 2)/(+, -, 0, n - 1) + (2n + 2) 2 /(+, - n - 1)) + 

/(+, 0, rc)(2n) 2 /(+, - n - 1) + /(+,-, n)(-2n) 2 /(+, +, n - 1) = 

t|^(3 + 2n) 2 (5 + 3n)(5 + 6n)(2n + l)! 2 . 

Then, taking a limit £/2 — > 0, wc obtain the contribution from this process to be 

28t 2 n(5 + 2n) 

~ 75{Uo + m) ■ 

Finally, 

28t 2 n(5 + 2n) 4(15 + 20n + 8n 2 ) 



£2 



75(J7 + 4£/ 2 ) 15(J7 + C/ 2 ) 



31 



APPENDIX B: DERIVATION OF THE EFFECTIVE MAGNETIC HAMILTONIAN FOR THE 

INSULATING STATE WITH TWO ATOMS 



To derive the effective Hamiltonian, we should be able to calculate matrix elements in J2TJ. Since energy and 
matrix elements in each subspace don't depend on z projection of total spin, we can choose S z components at our 
convenience. We can express any state \E\ >, \E$ > using known Clebsch-Gordon coefficients. For the state from 
\Eg > with S z = , we have 



= -2... 



\E,.S, = > C^^JN! = 2,Si = 2, S lz = m> \N 2 = 2, 5 2 = 2, S 2z = -m > 

.2 



Now we can write any one of the states \Ni, Si, Si Z > via creation and annihilation operators since we know how to 
express spin operators via creation and annihilation operators Ijlljl . Evaluation of e~6 — eg is quite simple since total 
spin conservation of the tunnelling term doesn't allow mixing of this subspaces with any other. Therefore, as in l|A2() . 
we just need to calculate the normalization of the state into which we hop. Using this procedure, we obtain energies 
(HU. 



Now let's consider energy in the \E\ >, \E 2 > subspace. From \Ei > we can hop only into high energy states 



|iVl : 
\Nl 



3, N 2 = 1, Si = S 2 = 1, Si + S 2 = >, 
= 1, N 2 = 3, 51 = 52 = 1,51 + 52 = >, 



(Bl) 



since in the Hilbert space of each well spin can change only by ±1. For N — 2 from \E 2 >, we can also tunnel only to 
these states since 3 and 1 cannot add to form total spin 0, and there is no state 5i = 52 = 3 which can also add up to 
total spin 0. Therefore, our exact Hamiltonian in the basis of \Ei >, \E 2 > and high energy states i|Bl|) has the form 









Vi 


Vi 





6LT 2 


v 2 


v 2 


Vi 


v 2 







\\ 


V-2 





U 



We can diagonalize this hamiltonian in the low energy \Ei >, \ E 2 > subspace in the limit 

Vi,V 2 < U ,U 2 < IV 

First, we integrate out high energy levels - this is done as described in [lj. We use the following matrix identity: 



= [{A-BD- 1 C)- 1 ] l3 , 



where 1 < i,j < 2. 

In our case, D has the form UqI 2 , so it is easy to calculate inverse matrix. Finally, our effective Hamiltonian has 
the form 







2 




ViV 2 


o m 2 




ViV 2 





Now we can diagonalize this 2x2 matrix; its energy levels are 

vi / 



3U 2 



V4 



U 



± \ (3U 3 



V 2 4- V 2 
Vl +V2 ) 2 + 12 



U 2 Vj 
U 



Using expressions for all states of interest in Fock basis, we can calculate 

^10 ~ 



Vi = -n 



3< y2 = -*V3' 



which leads to QQJ-^. 



32 



Now, let's calculate energy for the Si + Sj — 2 subspace. In this case we can hop to 4 states: 

|JVi = 3, N 2 = 1, Si = S 2 = 1, Si + S 2 = 2 >, 

|JVi = 1, JV a = 3, Si = S 2 = l,Si+S 2 = 2 >, 

\Ni = 1,N 2 = 3, Si = 1, S 2 = 3, S x + S 2 = 2 >, 



|JVi = 3, N 2 = 1, S 1 = 3,S 2 = l,Si+S 2 = 2> 
Matrix A has the form(in E 3 , E 4 , E§ basis) 

" w 2 

3U 2 
6C/ 2 



and matrix B has the form 



Vi v 2 v 3 
v 2 vi y 3 
v 4 v 5 y 5 



and C = -B T . The effective Hamiltonian is 



3U 2 
3Z7 2 
6C/ 2 



1 



2FiF 2 (Vi + V 2 )V 4 + V 3 V 5 

2v 1 v 2 vi + vi + vi {Vi + v 2 )v 4 + v 3 v 5 

(Vi + v 2 )v 4 + v 3 v 5 (v 1 + v 2 )v 4 + v 3 v 5 2(v 4 2 + vi) 



^From the symmetry of this matrix, one eigenvalue is easy to determine - it corresponds to eigenvector (1,— 1,0) T 
and equals 



3C/ 2 - 



Vj + Vj + Vj - 2V X V 2 



Two other eigenvalues are also easy to obtain -they are the solutions of some quadratic equation. Calculating Vi 



gives 



Energy levels are 



/2 /To /14 /14 12 



9U 2 - I2t 2 ± * / 144-^- + A0^-U 2 + 9U. 



U 



APPENDIX C: MEAN FIELD SOLUTION FOR THE CASE OF TWO BOSONS PER SITE 

Our mean field state depends on 12 variables - 6 complex numbers, subject to normalization H40(l . However, energy 
is clearly the same for all states that can be transformed into each other by global SU (2) rotations, so that gives us 
3 conditions we can choose. We also have an overall U{1) phase freedom, so the number of independent parameters 
reduces to 12 — 3 — 1 — 1 = 7. We can parameterize 



a —b e i<Pn 
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It is convenient to choose 3 constraints to be: 

b = 0,(p-i = tpi. 

We can always find an axis of quantization for which the first condition is satisfied and then we can satisfy the last 
condition on phases by rotating along chosen axis of quantization. We have 7 degrees of freedom: eight variables 
b± 2 , b±i, <p±2, fi, 9 with the normalization constraint 

b\ + b 2 _ 2 + b\ + b\ = 1. 

We can express mean field energy via these degrees of freedom, and it will be some polynomial up to the fourth order 
in 6 m (we won't explicitly write down this polynomial). First, let's make a transformation that diagonalizes the part 
quadratic in &± 2 , b±\. That orthogonal transformation is 

-u 2 + v 2 iti - v t 
6 - 2 = 2 ^ 



ui+vi , u 2 +v 2 

b 1 = ^^,b 2 = ^—. 

Now we solve the normalization constraint by 

vi = sin tp sin ipi,v 2 = cos ip sin ip 2 , 



u\ = sin ip cos tpi,U2 — cos tp cos ip 2 ■ 

The numerical procedure now consists of fixing the values of t and and minimizing the expression for energy over 
six angles using steepest descents method. Using this procedure we can numerically find energy as a function of 9 
for fixed t. If the minimum is attained at 9 = 0, then mean field wave function is a trivial singlet. Result of this 
minimization leads to the state 

1 1 1 1 T 

m ~ ( 76' V3' ' T 73' V6 } 

that can be rewritten as (|42l) after rotation. Though this result was obtained numerically, we can check analytically 
that at this point all first derivatives of the energy over angles ip, ip\ and ip2 vanish. 



APPENDIX D: LARGE TV EXPANSION 



In this Appendix we prove some properties of wave functions \ip) n <|<TOJl - The normalization factor J\f in <|fi9[) may 
be found by considering the overlap of two states. It is sufficient to consider wave functions constructed of only two 
single particle states since rotation in the spinor state can always bring our "pure condensate" wave functions to this 
case: 

\N,m) = i(cos(M++sMiaj,)*|o}, 

\N,n 2 ) = i(cos^2O++sin^2at) w |0). 

(Dl) 

We have 

(N, ni| TV, n 2 ) = (0| (cos 4>ia x + sin (f>\a v ) N {cos foa^ + shK/> 2 ap w |0) = 



1 

^E^) 2 ^ 08 ^ 008 ^)"^^! sin0 2 ) Ar - fc <O|a^- fc (4) fc (4) Ar - fc |O) = 

fe=0 
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-^^C , Ar( c °s0iCOS(/)2) fc (sin^isiii0 2 ) Ar k = ^cos N (fa - fa) = ^(n 1 n 2 ) N . 

fc=0 

Orthogonality and normalization for large N now become obvious after noting that for (111112) = cos9 and 8 < ir/2 
we have cos^ 9 « e~ Ne / 2 . 

To prove l|65[) we consider wave functions 

\N,n) = j^{n x ai+n y a\) N \Q), 

|JV + l,n') = -J—^ot+n^lO). (D2) 

Simple calculation gives 

(N + 1, n'\al\N, n) = ( f + 1)1 < (nnf = (JV + l) 1 / 2 n x M" - n'). (D3) 
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